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Although  the  study  of  wave  scattering  by  randomly  rough  surfaces  started 
several  decades  ago,  it  is  still  a  challenging  research  topic  today  [1-11]  . 
Basically,  when  the  surface  is  slightly  rough  and  its  surface  slope  is  generally 
smaller  than  unity,  the  perturbation  technique  can  be  used.  When  the  radius  of 
curvature  of  the  surface  is  much  greater  than  the  wavelength,  the  Kirchhoff 
approximation  can  be  applied  [12].  In  other  words,  for  relatively  long  wavelength 
waves,  the  perturbation  technique  is  a  good  choice;  for  relative  short  wavelength 
waves,  the  Kirchhoff  approximation  can  be  adopted.  In  any  situation,  second 
moments  of  the  scattered  wave  field  can  provide  some  useful  information  such  as 
average  intensity.  More  fruitful  information  on  scattering  characteristics  relies 
upon  the  knowledge  of  fourth  moments  of  the  scattered  wave  field.  However,  more 
difficulties  will  be  encountered  in  deriving  the  fourth  moments.  Wave  scattering 
by  randomly  rough  surfaces  is  not  only  an  interesting  topic  for  theoretical  study 
but  also  a  practical  problem  in  many  applications  ranging  from  microwave 
scattering  by  ground  surfaces  in  remote  sensing,  acoustic  scattering  by  the  ocean 
surface  and  floor  in  underwater  acoustics,  to  optical  scattering  by  a  rough  metal, 
dielectrics,  or  semiconductor  surfaces  in  designing  optics  and  electronics 
devices . 

In  this  paper,  three  types  of  fourth  moments  are  evaluated  for  acoustic  waves 
forward  scattered  by  a  rough  ocean  surface.  The  first  one  is  the  scintillation 
index  a  defined  by 


a  =  [<  |p|4>/<  |p|2>2  -1]1/2  . 


(1) 
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Here  p  is  the  scattered  acoustic  wave  field  and  <  >  stands  for  an  ensemble 
average.  The  scintillation  index  describes  the  intensity  fluctuation  of  the 
scattered  wave .  The  purpose  of  evaluating  the  other  two  fourth  moments  is  to 
extract  some  information  on  the  phase  fluctuation  of  the  scattered  wave.  To  see 
this  point,  let  us  first  consider  the  case  of  weak  scattering  in  which  the 
fluctuations  of  the  log-amp]  •'.tude  and  phase  of  the  scattered  wave  field  are  either 
weak  or  jointly  Gaussian  distributed.  Hence,  by  writing  the  scattered  wave  p  as 

(2) 

p  =  exp  (%  +  lS) 


in  terms  of  the  log-amplitude  x  and  phase  S,  the  two-position  coherence  function 
<(PiP2*)2>  is  given  by 

<(PiP2’)2>  =  exp  [2<  (Xi  +  Xa)2>]  exP  t~2<  ^S1-S2)2>] 
x  exp  [4i<  (Xi+Xj)  (■51-S2)  >]  , 


where  the  subscripts  1  and  2  stand  for  the  waves  at  positions  1  and  2, 
respectively.  Meanwhile,  the  two-position  intensity  correlation  function  <IX  I2> 
is  given  by 

<W  -  <!pi!2Ip,I2>  (< 

■  exp  [2 <  (Xi+X2>2>} 


where  the  intensity  I  -  | p | 2 .  Therefore, 

|<(PiPa)2>l 


<JiI2> 


=  exp  [  -  2  <  ( 5X  -  S2 ) 2  >  ] 


(5) 


which  describes  the  fluctuation  of  the  relative  phase  between  the  two  observation 
positions.  When  scattering  becomes  stronger,  the  simple  result  in  Eq .  (5)  is  not 
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valid.  However,  the  ratio  |<(p1p2*>|/<lil2>  or  the  difference  between  |<(PiP2*)2| 
and  <Iil2>  still  carries  information  on  the  relative  phase  fluctuations.  To 
extract  this  information,  the  normalized  two-position  coherence  function 
<(plp2*)2/<Ii><l2>  and  normalized  two-position  intensity  correlation  function 
<I1I2>/<I1>  <I2>  are  to  be  evaluated  in  this  paper. 

Because  the  acoustic  frequency  of  concern  is  relatively  high,  the  Kirchhoff 
approximation  will  be  used.  Two  types  of  spectral  functions  for  the  surface- 
height  fluctuation  are  considered:  a  Gaussian  spectrum  and  the  Donelan/Pierson 
(D/P)  spectrum.  The  latter  is  obtained  from  a  model  describing  the  fluctuations 
of  the  ocean  surface  height  which  are  controlled  by  the  wind  speed  on  the  ocean 
surface  [13].  Both  the  scale-size  and  mean-square  fluctuation  of  the  ocean 
surface  height  are  determined  by  the  wind  speed.  The  fluctuation  of  surface 
height  f  is  assumed  to  be  statistically  Gaussian  distribu  .sd.  Numerical 
techniques  will  be  designed  to  evaluate  the  multi -fold  integrations  of  the  fourth 
moments.  For  mathematical  simplicity,  only  two-dimensional  propagation  and 
scattering  are  considered.  The  rest  of  this  paper  is  organized  as  follows.  The 
geometry  of  problem  is  discussed  in  section  II.  Also,  the  average  intensity  of 
the  scattered  wave  field  is  evaluated.  The  derivations  and  numerical  results  for 
the  scintillation  index  are  given  in  section  III.  Sections  IV  and  V  are  devoted 
to  evaluation  of  the  two-position  coherence  function  and  two-position  correlation 
function,  respectively.  Conclusions  are  drawn  in  section  VI. 

II.  GEOMETRY  OF  PROBLEM 

The  basic  geometry  of  two-dimensional  scattering  by  a  randomly  scattered 
rough  ocean  surface  is  shown  in  Figure  la.  The  depth  and  smooth  ocean  surface  are 
defined  to  be  along  the  z  and  x  axes,  respectively.  A  line  source  and  a  line 
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receiver  are  located  at  certain  depths  so  that  the  incident  range  RA,  scattering 
range  R, ,  incident  angle  8it  and  scattering  angle  8 ,,  are  as  shown  in  the  figure. 
The  surface  height  fluctuation  from  the  mean  at  z  -  0  is  described  by  the  random 
field  J(x)  with  <{"(x)>  -  0.  By  ignoring  some  unimportant  factors,  the  wave  field 
at  the  receiver  in  the  Kirchhoff  approximation  can  be  expressed  by  [14] 


p  =  J  dx  exp  [-i  (x2/x2  +  2ox  +  2y{  ( x ) )  ] 


(6) 


where 

,  .  2  ,  »in2«i  .  sin^e.  (7) 

*  k  Rt  Ra 

a  =  -~  (cosdj  -  cosQs )  (8) 

and 

Y  =  +  sin0s)  .  (9) 

Cm 

Here  k  is  the  wave  number  which  is  assumed  to  be  a  constant  in  the  ocean.  In 
other  words,  the  ocean  is  assumed  to  be  a  homogeneous  medium  to  the  sound  wave. 

In  Eq.  (6),  since  f(x)  is  a  random  field,  the  wave  field  p  is  stochastic.  Note 
that  the  coordinate  center  can  be  properly  chosen  so  that  8t  -  8t  «■  V  and,  hence, 
a  -  0. 

The  average  intensity  of  the  scattered  wave  field  can  be  easily  computed  to 

give 


<  I>  =  <  jp[2>  =  Ttxf 


(10) 
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Hence,  the  average  intensity  has  nothing  to  do  with  random  fluctuations  of  surface 
height. 

Figure  lb  shows  the  geometry  for  the  evaluation  of  the  two-position  intensity 
correlation  function.  Two  line  receivers  are  placed  at  different  depths  but  the 
same  horizontal  positions.  The  depth  difference  or  separation  is  d.  The  ranges 
Rx  and  R2  and  the  corresponding  scattering  angles  6i  and  $2  are  defined  in  the 
figure.  The  range  R,  and  the  corresponding  angle  6  are  used  to  define  the 
midpoint  between  the  two  receivers.  Hence,  Rlt  R2,  8  ± ,  and  82  can  be  expressed  in 
terms  of  R,,  8,  and  d  as 

0X  =  0  -  sin'Md  cos0/ 2R± ) 

02  =  0  +  sin'Md  cos0/2fl2) 

and 

ri.2  =  +  (d/2)2  t  d  Ra sin0]1/2 


(12) 

(13) 


These  equations  will  be  used  in  the  following  sections. 

III.  SCINTILLATION  INDEX 

Since  the  scintillation  index  is  related  to  a  limiting  value  of  the  two- 
positfon  intensity  correlation  function  <I1I2>,  the  derivation  in  this  section 
starts  with  <1  il2> -  Only  the  detailed  results  of  the  scintillation  index  are 
presented  in  this  section.  Those  of  the  two-position  intensity  correlation 
function  will  be  given  in  t.e  next  section. 

The  derivation  starts  with 
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<i1j2>  -  <|p1|j|p2|2> 

•• 

=ffff  dx1dxidx2dxA  exp {-i  [  (x2  -  x\ )/xjl  +  (x2  -  x42)  / x\2 

-to 

+  2a1  (xx  -  x2)  ♦  2a2  (x3  -  x4)]} 
x  <exp(-i[2Ya(((xi)  -  {(x2)  )  ♦  2y2(C(*3)  ~((x4)]}> 


The  ensemble  average  can  be  reduced  to 


exp  ( -4<(2>H) 

=  exp  {-4<C2>  [Yi  +  vl  ~  Yi (x:  -  Xj)  -  Y2Cc(x3  -  x4) 

+Y1Y2  (C{(x 1  -  x3)  -  C{(x2  -  x3)  -  C( (x2  -  x4)  +  Cc(x2  -  x4))]) 


Here,  <f2>  is  the  mean  square  fluctuation  of  the  surface  height  and  Cr(x)  is  the 
normalized  correlation  function  of  the  surface  height  fluctuation.  In  Eqs .  (14) 
and  (15),  xn ,  xf2,  7lt  and  y2  are  defined  in  Eqs.  (7)  and  (9)  for  the  two 
receiving  points,  respectively.  By  considering  the  following  variable 
transformation  with  unit  Jacobian: 


*1 

llll' 

t 

*2 

1-1  1-1 

P/2 

*3 

1-1-1  1 

S/2 

*4 

1  1  -1  -1 

.ff/4 

(16) 


the  integrations  with  respect  to  t  and  q  can  be  easily  completed.  The  result  is 


<I1  I2>  ff  c?Pcf5  exp  {-i  [2P~S/  (Ax  +  A2)  +  PB1] } 


(17) 


X  exp  [-4<C2>H(5,P)] 
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with 

H(S,P)  =  vl  +  Ya  -  Y?Cc(2PCri)  -  ylcz  (2*J(1  -  ii) ) 

+  YiY2fc«<^  +  3)f)  -  C<(ft(2n  -  1)  -31)  ~  C((P«(2n  -  1)  +  SC)  <18> 
+  Cc(CS-P)  0] 

in  terras  of 


Ali2  =  xfi.  2/2  C2 


(19) 


fl;  =  4C(a1x|J  -  a2Xf2)  /  (x|j  ♦  xj2) 


(20) 


T1  =  Xfi/  (Xfi  ♦  Xf2) 


(21) 


the  3  «  P/«,  3  =  S/(.  The  notation  i  represents  the  characteristic  length  of 
surface  height  fluctuation  or  more  precisely  is  the  scale  size  of  the  correlation 
function  Cf(x).  If  we  define 

A  =  X//2C2  <22> 

and  xf  and  7  as  in  Eqs.  (7)  and  (9)  for  the  midpoint  between  the  two  receivers, 
<I1I2>  in  Eqs.  (17)  and  (18)  can  be  further  reduced  to 


<Va> 


4n  C4AA1A2 
A2  +  A2 


J/  dPf  dS' exp[-i(DP'S'  +  BP')} 


x  exp  [-QZH'  (S',  P ') ) 


(23) 


with 
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H'lS'.P')  ~  Yi  ^  Yi  “  YiC  (2P7T))  -  y^C(2P/ (1  -T))  ) 
+  YiY2  IC(P'  +  S')  -C(P/(2»i-1)  -  S') 

-c(p7  (2n  - 1)  +st)  +  c(s/  -  p1)  ] 


Here, 


D  =  2A/  <At  +  Aj) 

(25) 

B  =  y/KB' 

(26) 

Yi,2  =  Yi.2/Y 

(27) 

4>2  =  4<C2>y2 

(28) 

C(x)  =  Cz(xx(/yf2) 

(29) 

Also,  the  new  integration  variables  are  P1  •  P/  (A) J/a  and  S1  -  3/  (A) 1/2 .  From  Eq. 
(28),  4>  is  an  indicator  of  scattering  strength. 

Equations  (23)  and  (24)  will  be  used  for  evaluating  the  two-position 
intensity  correlation  function.  The  rest  of  this  section  is  devoted  to  the 
computations  of  the  scintillation  index.  When  positions  1  and  2  coincide, 

B  -  0,  rj  -  1/2,  and  D  -  1,  as  can  be  seen  in  Eqs.  (20),  (21),  and  (25), 
respectively.  In  this  situation,  Eqs.  (23)  and  (24)  become 

m 

<I2>  ~  2nt*A2 J J  dP1  cfs'exp  ( -i  P'S')  exp  [ -<b2H'  (S',  p') )  <3°) 


with 


H'iS'.P')  =  2-2 C(P')  -2 C(S')  +  C(P/  +  5/)  +C(P/-5/) 


(31) 
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The  fact  that  the  normalized  correlation  function  Cf(x)  or  C(x)  is  an  even 
function  has  been  used.  It  is  noted  that  the  integration  in  Eq.  (30)  must  result 
in  a  real  number  for  <I2>  although  the  integrand  is  complex.  Hence,  this 
integration  possesses  inversion  symmetry  with  respect  to  both  P'  and  S'. 

Therefore,  Eq.  (30)  can  be  rewritten  as 

m 

<I2>/<±,2  =  —  Re  ff  dp'ds'exp  (-ip's1)  exp  [-<&2tf'  (S',  P1)  ]  (32> 

71  o 

where  <I>  given  in  Eq.  (10)  is  used  and  Re  stands  for  "the  real  part  of." 

To  numerically  evaluate  the  integration  in  Eq .  (32),  the  integration  area  is 
divided  into  three  regions  (I,  II,  III)  as  shown  in  Figure  2.  The  dimensions 
defining  the  three  parts  L,  Lj,  L2  and  L3,  follow  the  relations 

L2  *  +  sflL  <33> 

and 

Lj  =  yflLx  +  L  <34> 

The  sizes  of  L  and  Lx  will  be  defined  later.  In  these  three  integration  regions, 
different  approximations  will  be  used  for  numerical  computations.  If  the  physical 
size  of  L  is  much  larger  than  the  scale  size  of  the  correlation  of  f,  i.e., 

C(L)  «1 ,  the  integration  for  <I2>/<I>2  over  region  I  becomes 

(  <!2>/<I>2)  ,  =  —■  Re  f  f  dP'dS'exp  (-iP's')  exp  [-<b2H'  (S1,  P') }  (35) 

n  J  J  (i) 


with 


H1  {S',  P1)  ■  2  +  C(P'  -  S')  . 


(36) 
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If  the  variable  transformations  x  -  (S'  +  P')/2  and  y 
easy  to  obtain 

P'S'  =  x 2  -y2/ 4. 


P'  -  S'  are  used,  it  is 


(37) 


With  the  phase  factor  y2/4  in  the  integrand  of  Eq.  (35),  the  contribution  to  the 
integration  is  negligible  for  large  y.  Hence,  if  Li  is  chosen  so  that  L12/4»2rr, 
Eq.  (35)  can  be  approximated  by 


( <J2>/<I>2)  r  «  —  Re  j  dxj  dyexp[-i  (x2 -y2/4)  ] 

•Lj/y/'T 


x  exp  [-02  (2  +  C(y) )  ]  , 


(38) 


Next,  if  #  is  not  extremely  large,  we  expect  that  1  -  exp  [ -t2  C(y)]  decreases 
very  fast  with  y.  Therefore, 

M 

(<j2>/<I>2)r  =  ■Z-ReZ2(2ll  -  f  dy[  1  -  exp (-$2C(y) )  ]} 

-m 

x  exp  ( -2®2) 


where 


(1  +  i) 


(AO) 


and  £2  is  a  Fresnel  integral  as 

m 

S2  =  j  dxexp  ( -ix2)  . 

Ll/Si 


(41) 


The  integrations  in  Eqs.  (39)  and  (41)  can  be  easily  completed  for 
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<  <I2>/<I>2)r. 

For  region  II,  we  have 

•  L 

(<l2>/<!>2)  IZ  =  —  RejdP'jdS'exp  (-is'P')  exp  [~Q2h'  (s' ,  P1) }  .  <*2> 

71  L,  0 

Here,  a  factor  of  2  has  been  Included  since  region  III  covers  two  separate  parts 
which  result  In  the  same  Integration  value.  Because  the  physical  size  of  L  Is 
much  larger  than  che  scale  size  of  C(x),  the  upper  integration  limit  L  for  S'  in 
Eq .  (42)  can  be  replaced  by  <»  without  significantly  changing  the  integration 
result.  Also,  H'(S',P')  can  be  approximated  by 

H'  (S1,  P')  =  2  -  2C{S')  .  (43) 


Hence,  Eq.  (42)  becomes 

••  m 

(<I2>/<I>2)  ZI  =  —  Re  f  dP' f  dS'exp  (-is'P1)  exp  [-2®2(1  -  C(S') )  ]  .  («*) 

It  is  evident  that  the  integration  with  respect  to  S '  is  a  Fourier  transfon  and 
can  be  easily  evaluated  using  the  Fast  Fourier  Transform  algorithm  on  a  computer. 

For  region  III,  no  approximation  can  be  made.  Double  integrations  on  a 
computer  must  be  performed.  The  contribution  from  region  III  is 

L,  L  *  Lj  •  p1  i  l*  l,  -  p' 

(<I2>/<I>2)  zxz  =  —  J?e[ f  dP'  f  ds'-fdp'  f  dS '} 

0  0  0  Lj  '  ' 

x  exp  ( -is1  P')  exp  [~Q2H'  (S',  P')  ] 
with  H' (S' ,P' )  given  in  Eq.  (31). 

Hence,  <I2>/<I>2  is  the  summation  of  the  results  in  Eqs .  (39),  (44),  and 
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(45).  For  numerical  computations,  the  correlation  function  of  surface  height 
fluctuation  f  must  be  chosen.  To  illustrate  the  dependence  of  the  scintillation 
index  on  the  parameters  ♦  and  A,  a  Gaussian  correlation  function,  corresponding  to 
a  Gaussian  spectrum,  is  used 

C((x)  =  exp  [-  (x/t)2]  <46> 

and,  hence 

C(x')  =  exp  [-  (x'A)2]  •  (47) 


In  numerical  computations,  L  -  5 />/R  and  L i  =  (80n)l/2  are  used.  The  numerical 
accuracy  was  checked  by  decreasing  the  step  sizes  in  the  integrations  until  the 
results  were  not  changed,  Three  curves  are  plotted  in  Figure  3  for  the 
scintillation  index  a  versus  $  for  three  different  A  values  of  2.58  (long-dashed), 
50  (short-dashed),  and  245  (solid)  when  the  Gaussian  correlation  given  in 
Eqs .  (46)  and  (47)  is  used.  Each  curve  increases  with  4>  and  approaches  unity 
asymptotically.  In  other  words,  the  scintillation  index  increases  with  scattering 
strength  and  becomes  close  to  one  when  saturation  is  almost' reached.  The  wave 
field  follows  a  jointly  Gaussian  distribution  in  the  saturation  regime.  The  fact 
that  the  scintillation  index  does  not  exceed  unity  implies  that  the  phenomenon  of 
focusing-defocusing  does  not  occur.  Among  these  three  curves,  the  higher  the 
value  of  A,  the  larger  is  the  scintillation  index.  To  further  explore  this 
trend,  from  Eqs.  (7)  and  (22)  under  the  assumptions  8t  ~  8t  —  6  and  Rt  —  R,  -  R,  we 
can  obtain 

\  =  R/ (2 kt2 s in26)  .  (A8) 
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It  is  easy  to  understand  that  a  smaller  i  value  (larger  A)  leads  to  stronger 
scattering  and,  hence,  a  larger  scintillation  index.  Also,  as  the  range  R  is 
increased  more  complete  evolution  from  phase  fluctuations,  which  are  due  to  rough 
surface  scattering,  into  amplitude  fluctuations  is  expected.  Therefore,  a  larger 
R  value  (larger  A)  results  in  a  larger  scintillation  index.  It  is  noted  that  for 
<t<2 ,  the  approximations  used  fail.  The  numerical  methods  used  will  also  fail  for 
very  strong  scatter. 


IV.  TWO-POSITION  CORRELATION  FUNCTION 

The  computations  for  the  two-position  correlation  function  start  from 
Eqs .  (23)  and  (24).  By  using  Eq.  (10), 

«•  «t 

<  J,J2> /<!,><  J,>  =  —  —  ft-:---  Re  f  dP1  [  dS' exp [ -i  (DP1  S'  +  BP1)  ] 

n  (At  +  A2)  JQ 

X  exp  [-Q2H'  (S',  P1)] 


where  H'(S',P')  is  given  in  Eq.  (24).  In  obtaining  Eq .  (49),  the  symmetry  with 
respect  to  P'  was  used.  Because  this  integration  does  not  possess  symmetry  with 
respect  to  S',  the  partition  of  the  integration  area  is  different  from  that  for 
computing  the  scintillation  index.  For  the  integrations  in  Eq .  (49),  four  regions 
(I  through  IV)  are  designated  as  shown  in  Figure  4.  Among  them,  regions  I  and  II 
individually  have  two  separate  sections.  For  simplicity  in  notation,  we  define 


E  =  2 Re  f  dP' I  dS'  exp  [-i  {DP1  S'  +  BP') ) 


x  exp  [-$ 2H,{S,I  P ')  ] 


and,  hence 
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<I1I2>/<I1XIJ> 


A E 

u  (Ax  +  A2) 


(51) 


The  dimensions  L,  Llt  L2,  and  Lj  in  Figure  4  are  the  same  as  those  in  Figure  2. 

The  approximation  used  for  region  1  is  the  same  as  that  for  region  I  in 
Figure  2.  It  results  in 

Er  =  2i?e(F‘G)  (52) 

where 


F  =  <  — ),/2exp  (-iB2/4D){[l -CAu)  -(Mv)]  -  i  [1  -SAu)  -5i(v)]}  (53) 

2D 

and 


G  =  exp  [-G2  (y2  +  y?)  ]  {( ) 1/2  (1  +  i)  exp  (iflJ/4I>) 

-2  J  dy[  1  -  exp  (-^YiYjCfy) )  ]} . 

o 


(54) 


The  approximation  for  region  II  here  is  again  the  same  as  that  for  region  II  in 
Figure  2.  The  contribution  from  this  region  is 


En  =  Re  f  dS' J  dP1  exp  ( -iP'S') 

DLj  — 

x  (cos  (P'B)  [exp{-<tJ[y|  +y2  -  y2C(2P/T])  -y^C(2P/(l  -r\)  )  ]  ) 
-exp  ( -<J>2  ffj  +  y2) )  3). 


(55) 


Again,  the  algorithm  of  the  Fast  Fourier  Transform  can  be  used  for  the  integration 
with  respect  to  P'.  No  approximation  can  be  made  for  region  III;  its  contribution 
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is 


L,  L  +  L>-P'  i  L, 

EIU  -  4Re[f  dP/  f  dS'  +  fdPffdS']  exp  (-iDP'S1) 

L  0  0  0 

X  cos  ( P'B )  exp  [-<1 >2H'(S',  P')  ] 


(56) 


Here,  H'(S'.P')  was  defined  by  Eq.  (24).  Finally,  for  region  IV  the  two-fold 
integration  cannot  be  simplified.  However,  H'(S',P')  can  be  reduced  to  produce 


EIV  =  2  Ref  dP'  (L*''P'~L')  <2n'1)  dS'  exp  (-iDP'S'  -  iP'B) 

J  J-l-ip'-Lj)  (2n-i> 

x  exp{-$2  [Yi  *  71  “  Y1Y2  ( C(p '  ( 2t|  -  1)  -S')  *  CiP1  ( 2rj  -  1)  +  S')  )  ] ) 


This  two -fold  integration  is  quite  time  consuming  on  a  computer,  especially  when  rj 
approaches  1/2.  It  is  noted  that  E!V  is  equal  to  En  when  *7  -  1/2.  This  equality 
can  be  used  to  check  the  accuracy  of  EjV.  Combining  Eqs.  (52),  (55),  (56),  and 
(57),  we  can  obtain  E  as 

E  *  El  +  E11  +  EUI  -  E ^  (58) 


and  thence  <IlI2>/<I1Xl2>  from  Eq.  (51). 

To  describe  realistic  ocean  surface  fluctuations,  the  D/P  spectrum  [13)  is 
used  at  low  wave  numbers  and  the  empirical  spectral  model  of  Pierson  and  Stacy 
[15]  and  Pierson  [16]  at  high  wave  numbers.  This  combined  spectrum  has  an 
approximate  dependence  on  the  inverse  cube  of  the  wave  number  over  a  wave  number 
span  that  is  determined  by  wind  speed  and  vanishes  outside  this  region.  In  this 
spectrum,  both  scale  size  i  and  mean-square  fluctuation  <f2>  are  controlled  by  the 
wind  speed  at  the  ocean  surface.  Two  normalized  D/P  correlation  functions  are 
depicted  in  Figure  5  for  wind  speeds  of  10  m/s  (solid  curve)  and  15  m/s  (dashed 
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curve).  Side  lobes  f.  •.  both  curves  can  be  seen.  The  scale  size  i  is  defined  as 

the  length  at  which  the  correlation  drops  to  one-half  of  the  maximum  value. 

Therefore,  we  have  i  -  9.27  m  and  22.61  m  for  wind  speeds  of  10  m/s  and  15  m/s, 
respectively.  Also,  the  mean- square  fluctuations  are  <f2>  -  1.08  m2  and  5.38  m2, 
respectively,  for  the  lower  and  higher  wind  speeds. 

In  numerical  computations,  we  again  choose  L  -  5/ >/K  and  Lx  -  (80n)1/2 

Other  parameters  are  9  -  10°  and  k  -  (frequency  -  3  kHz).  Two  values  for 
Rj  -  R,  -  R  at  1  km  at  10  b  will  be  used.  Figure  6  shows  the  results  for  the 
normalized  two-position  intensity  correlation  function  <I1I2>/<I1XI2>  as  a 
function  of  the  separation  d  between  the  two  receivers  for  /arloas  situations. 

The  four  solid  curves  labelled  by  A,  B,  C,  and  D  show  the  results  for  the  )/P 
correlations  with:  (A)  R  -  1  km  and  wind  speed  -  10  m/s, 

(B)  R  -  1  km  and  wind  speed  -  15  ra/s, 

(C)  R  -  10  km  and  wind  speed  -  10  m/s, 

and  (D)  R  -  10  km  and  wind  speed  -  15  m/s. 

For  comparison,  two  dashed  curves  are  plotted  for  the  Gaussian  correlation  with 
(A')  R  -  1  km  and  (C')  R  -  10  km.  The  scale  size  of  the  Gaussian  correlation  is 
set  at  9.27  m  and  mean-square  fluctuation  <f2>  is  1.08  m2  which  correspond  to  the 
wind  speed  at  10  m/s  in  the  D/P  correlation.  The  comparison  between  curves  A  and 
B  shows  that  for  weak  scattering  the  two-position  intensity  correlation  is  higher 
for  a  higher  wind  speed.  When  the  range  R  increases,  phase  fluctuations  of  the 
scattered  wave  can  evolve  into  amplitude  fluctuations  more  completely  and,  hence, 
the  curves  for  the  intensity  correlation  become  higher.  In  the  case  of  R  -  10  km, 
although  the  scintillation  index  at  a  wind  speed  of  15  m/s  is  higher  than  that  at 
10  m/s,  the  intensity  correlation  decreases  faster  with  separation  between  the  two 
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receivers.  The  comparison  between  curves  A  and  A'  and  that  between  curves  C  and 
C'  show  that  if  all  other  parameters  are  the  same,  the  Gaussian  spectrum  leads  to 
weaker  scattering.  Because  forward  scatter  in  the  Kirchhoff  approximation  is 
governed  by  local  specular  reflection  from  properly  oriented  "facets"  on  the 
surface,  the  mean-square  surface  slope  determines  the  area  of  the  surface  that 
contributes  effectively  to  scattering.  The  quasi  power  law  dependence  of  the  D/P 
spectrum  yields  higher  mean-square  surface  slopes,  and,  hence,  stronger  scattering 
than  does  the  Gaussian  spectrum.  It  is  noted  that  curves  A,  B,  C,  and  C'  are 
almost  parallel,  indicating  that  the  correlation  lengths  in  these  cases  are  about 
the  same.  Since  curve  D  is  steeper  than  curve  C,  the  correlation  length  is 
shorter  and  scattering  is  stronger  in  the  case  of  D.  For  the  same  reason,  the 
scattering  in  the  case  of  A’  is  very  weak. 

V.  TWO -POSITION  COHERENCE  FUNCTION 

As  discussed  earlier  in  this  paper,  the  two-position  coherence  function 
carries  information  on  phase  fluctuations.  The  derivations  start  with 

m 

<  (P1P2 )  2>  =  ////  dx1dx2dx3dxt  exp (-i  [  (x2  •*-  x32 )  /xjx 

(59) 

-  (x22  +  x2)  /Xf2  +  2a:  (x:  +  x3)  -  2a2  (x2  +  x4)  ] ) 
x  <exp{-i  [ 2 Yi  (C  (*i)  *  C  (x3) )  -  2y2  (C  (x2)  +  f  (x4)  )  ]  }> 

After  using  the  variable  transformation  in  Eq .  (16),  the  integrations  with  respect 
to  t  and  q  can  be  completed  using  the  method  of  stationary  phase  to  produce 
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<  (PiPa*>  2> .  -  2n  xflx(3fexp  [~2±  (a.\x}2  -  afx^)  ] 

x  fdPf  cf5exp(-i  [  ("P +  3)2/4A1  -  (P-3)2/4Aj]> 

O  -• 

x  exp  [-4<{2>//(P,P)  ] 


where 

tf(S,P)  =  Yi  +  Yi  +  Y?Cc(P  +  3)  +  Ya  C{  (3  -  P) 

-YxY,IC{  (?*^/2)  +C(CS  +  T0/2)  +CcCq0/2-S) 

+cc(T0/2  -p)l 


with 


q0  =  4«2  (a2A2  -  a^) 


(62) 


Further  normalization  of  P  and  3  does  not  simplify  the  computations.  Again,  for 
simplicity  of  notation  the  two-fold  integration  A  is  defined  by 


A  =  fdPf  d5  exp(-i  [  (P  +  3)  2/4A1  -  (P -3)  2/4A*] } 

O  -«• 

x  exp  [-4<C2> HCs.’P)  ] 


(63) 


Hence,  the  normalized  two-position  coherence  function  is 

<  (p1Pj)2>/<I1XI2>  =  — -  exp  [~2i  [dlxU  -  tt?x?,)  1  (64) 

Xy/AiAj 

Five  regions  can  be  identified  for  computing  A  as  shown  in  Figure  7.  Regions 
II  and  V  individually  include  two  areas.  The  dimensions,  L' ,  ,  L2' ,  and  L3' 

are  clearly  shown  in  the  figure.  They  are  defined  by 
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(65) 


and 


L3  =  v^L/  +  L'. 


(66) 


The  choice  of  L'  and  Lj '  must  satisfy  the  inequalities  L1  -  g0/2>  >1  and 
L'  1*/4AJ> >2n  for  both  j  -  1  and  2,  respectively.  For  the  contributions  from 
regions  I  and  IV,  the  approximation  used  for  region  I  in  Figure  2  can  be  applied. 
The  results  are 


Ax  =  exp  [-®2  (yj  +  y!>  1  9^  (2*^) 1/2  (1  +  i) 

-2  J dy[l  -  exp  (-$2y?Cc(y) )  ]} 

O 


(67) 


and 


AIV  =  exp  [-#2(Yi  +Yz)]  g,{(2irA1) 1/2  (1  -i) 

m 

-2  j dy [1  -exp  (-4>2y! Cc(y) )  ]} 

O 


(68) 


where 


=  (nA1/2)1/2(  [1/2  -1(1/2  -5^1^//^)])  <69) 


and 


g2  =  (JtA2/2)1/2([l/2  -C^lU^^))  +  i[l/2 


(70) 


21 


9  August  1990 
CCY: STM: lzh 


To  obtain  the  contributions  from  regions  II  and  V,  the  approximation  used  for 
region  II  in  Figure  2  can  be  utilized  to  produce  the  combined  result 

aii  v  “  A xx  +  Av 


4AA 

Al+A2 


m 

/ 

(A*  *A,)/2AiAj 


dPex p  [-iAiAj  -  Ax)  P5/  (Ax  +  Aj) 2] 


x  J  dSexp  (-i?3)  exp  [-i  (Aj  -  Ax)  S7/ (dA^)  ] 


(71) 


X  (exp{-4>2  [Yi  +Y2 -YiYs  (C<  (3  +  <70/2)  +  C(  (5  -  q0/2)  )  ] } 
-  exp(-4>2(Yi  +Yz) ))  • 


Again,  the  FFT  algorithm  can  be  used  for  the  integration  with  respect  to  3. 
Finally,  without  any  approximation  the  contribution  from  region  III  is 

L'  +  L'  +  Li -7  l'.L,' 

AIII  =  [  /  cfP  /  dS  -2  j  dp  f  ds] 

°  Li  -(L'*Li-?)  (72) 

x  exp{-i  [  (P+5)2/4A1  -  (P-'3):>/4A2]}exp  [~<b2H'  (S,  P)  ] 

where  H1  (3,P)  is  the  same  as  P(3,P)  in  Eq.  (61)  except  that  7i  and  y2  are 

replaced  by  77  and  Ya  •  respectively.  The  normalized  two-position  coherence 
function  can  be  obtained  from  the  equality  A  -  A:  +  AIV  +  An  v  +  Ain  and  Eq .  (64). 

The  same  parameter  values  as  before  are  used  for  numerical  evaluation  of  the 
normalized  two-position  coherence  function.  L’  and  are  chosen  so  that 
L'  -  ~q~0/2  -  5  and  l[  =  (80nA)J/2.  Figure  8  shows  the  absolute  values  of  the 
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normalized  two-position  coherence  function  |<(PiP2*)2>|/<ii><i2>I  as  a  function  of 
the  separation  d  for  various  situations  which  are  the  same  as  those  in  Figure  6. 
Each  curve  drops  from  its  maximum  value,  which  is  <I2>/<I>2,  to  zero.  All  the 
letters  indicate  the  same  cases,  respectively,  as  those  in  Figure  6.  The  shorter 
coherence  lengths  in  the  cases  of  B  and  D  when  compared  with  the  cases  of  A  and  C 
confirm  that  the  scattering  at  wind  speed  15  m/s  is  stronger.  Apparently,  phase 
fluctuations  play  an  important  role  in  determining  the  size  of  coherence  length. 
Curves  A'  and  C'  show  that  the  coherence  lengths  in  the  case  of  the  Gaussian 
spectrum  are  about  the  same  as  those  in  the  case  of  D/P  spectrum  when  scattering 
is  weak  and  are  smaller  when  scattering  is  stronger. 

As  discussed  in  section  I,  when  scattering  is  weak,  the  ratio  between  the 
absolute  value  of  the  two -position  coherence  function  and  the  two-position 
intensity  correlation  describes  the  fluctuation  of  the  relative  phase  between  the 
two  positions  [see  Eq.  (5)).  Since  the  scintillation  index  is  0.52  in  the  case 
when  R  -  1  km  and  the  wind  speed  is  10  ra/s  in  the  D/P  spectrum,  Eq.  (5)  must  be 
approximately  true.  By  using  Eq .  (5),  the  root-mean-cquare  of  the  fluctuation  of 
relative  phase,  i.e.,  <(Sj  -  S2)2>1/z,  is  plotted  for  this  case  in  Figure  9.  The 
fluctuation  of  the  relative  phase  almost  reaches  2.25  (radian)  at  d  -  17.5  cm, 
which  is  only  0.35  times  a  wavelength. 

VI .  CONCLUSIONS 

Three  types  of  fourth  moments  of  forward-scattered  acoustic  waves  from  a 
randomly -rough  ocean  surface  have  been  evaluated.  The  first  is  the  scintillation 
index  which  characterizes  the  intensity  fluctuations  of  the  scattered  wave.  The 
second  is  the  two-position  intensity  correlation  function.  It  indicates  the 
spatial  correlation  of  the  intensity.  The  third  is  the  two-position  coherence 
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function  which  carries  information  on  phase  fluctuations  of  the  scattered  waves. 
Particularly,  when  the  scattering  is  not  strong,  the  ratio  of  the  absolute  value 
of  the  two-position  coherence  function  over  the  two-position  intensity  correlation 
Sanction  exactly  describes  the  mean- square  fluctuation  of  the  relative  phase 
between  two  observation  positions. 

The  Fresnel-corrected  Kirchhoff  approximation  was  used  to  obtain  integral 
expressions  for  the  fourth  moments.  Various  approximation  techniques  were 
developed  for  numerically  evaluating  the  integrals.  The  approximations  and 
numerical  methods  are  not  applicable,  however,  for  treating  the  very  weak  or 
extremely  strong  scattering  regimes.  Two  types  of  power  spectra  for  the  surface 
height  fluctuation  were  considered:  a  Gaussian  spectrum  and  the  Donelan/Pierson 
spectrum  which  is  an  empirical  model  based  on  ocean  wave  measurements.  The  D/P 
spectrum  resulted  in  stronger  scattering  than  the  Gaussian  spectrum  for  the  same 
values  of  A  and  *  which  was  attributed  to  the  greater  mean-square  surface  slope 
obtained  with  the  D/P  spectrum.  The  two-position  coherence  function  was  found  to 
decay  much  more  rapidly  with  displacement  between  receivers  than  the  intensity 
correlation  function.  For  weak  scattering,  this  led  to  large  relative  phase 
fluctuations  between  two  vertically  displaced  observation  points. 
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Figure  1(a).  Geometry  of  problem.  The  acoustic  wave  is  transmitted  from  a  line  source,  scattered  by  the 
rough  surface,  and  received  by  a  line  receiver.  No  variation  for  acoustic  waves  along  the 
y  direction  is  assumed.  Surface  height  fluctuation  is  denoted  by  Part  (a)  shows  a  general 

geometry . 


Geometry  of  problem.  The  acoustic  wave  is  transmitted  from  a  line  source,  scattered  by  the 
rough  surface,  and  received  by  a  line  receiver.  No  variation  for  acoustic  waves  along  the 
y  direction  is  assumed.  Surface  height  fluctuation  is  denoted  by  £.  Part  (b)  demonstrates 
explicitly  the  geometry  of  two  line  receivers  separated  by  d  in  depth. 


correlation  for  three  A  values: 


P'(2  7)  -1)=-S*  p  P'(2  T)  -1  )=S’ 
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Figure  Four  regions,  1  through  IV  in  the  first  and  second  quadrants  of  the  S'  -  P*  plane,  are  designated 

for  the  integrations  for  the  two— position  intensity  correlation  function.  Regions  I  and  II 
individually  include  two  separate  areas. 


31 


9  August  1990 
CCY:STM:lzh 


© 


(xpo  uoi.;B|ejjoo 


x  (m) 

Figure  5.  Correlation  functions  of  surface  height  fluctuation  C  (x)  in  the  Donelon/Pierson  model  when 
the  wind  speeds  on  the  ocean  surface  are  10  m/s  (solid  curve)  and  15  m/s  (dashed  curve). 
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Figure  9.  Root-mean-square  of  the  fluctuation  of  the  relative  phase  between  two  receivers  in  the  case 
R  =  1  km  for  a  wind  speed  of  10  m/s  and  the  D/P  spectrum. 
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